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We solve for single distorted black hole initial data using the puncture method, where the Hamil- 
tonian constraint is written as an elliptic equation in _R 3 for the nonsingular part of the metric 
conformal factor. With this approach we can generate isometric and non-isometric black hole data. 
For the isometric case, our data are directly comparable to those obtained by Bernstein et at, 
who impose isometry boundary conditions at the black hole throat. Our numerical simulations are 
performed using a parallel multigrid elliptic equation solver with adaptive mesh refinement. Mesh 
refinement allows us to use high resolution around the black hole while keeping the grid boundaries 
far away in the asymptotic region. 
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I. INTRODUCTION 

Distorted black holes are expected to be produced 
by astrophysical events such as asymmetrical gravita- 
tional collapse, black hole-black hole coalescence, black 
hole-neutron star coalescence, and possibly neutron star- 
neutron star coalescence. As ground based gravitational 
wave detectors such as LIGO, VIRGO, TAMA, and GEO 
increase their sensitivities, and the planned launch date 
for the space based detector LISA grows near, the need 
for researchers to develop a theoretical understanding of 
these systems becomes increasingly important. In sup- 
port of that effort, we need to develop the techniques and 
tools to evolve a dynamic, distorted black hole and ex- 
tract the emitted gravitational radiation as it settles to 
a quiescent state. 

One of the difficulties encountered in the numerical 
treatment of single and multiple black hole systems is 
the presence of nontrivial topologies. In the past, re- 
searchers have addressed this problem by using isometry 
conditions at inner boundaries to represent black hole 
throats. Another approach to black hole evolution is ex- 
cision, in which a region inside each apparent horizon 
is removed from the computational domain. Recently it 
has been shown that black holes can be treated in terms 
of fields on R 3 by splitting the conformal factor for the 
spatial metric into singular and non-singular terms. This 
so-called "puncture method" was applied to the Bowen- 
York pj family of black hole initial data sets by Brandt 
and Briigmann Q, and used for evolution studies by 
Briigmann 0] and others (see, for example, Refs. 0,0)- 

In this paper we make a modest extension of the punc- 
ture construction for initial data to include single dis- 
torted black holes. We reproduce and extend the results 
of Bernstein et al. 0, 0, Ig , who constructed "black hole 
plus Brill wave" initial data sets using isometry condi- 
tions at the black hole throat. Whereas the black holes 
obtained by Bernstein et al. are isometric by construc- 
tion, our distorted puncture black hole data sets are iso- 
metric only when the free parameters \x and to, defined in 
Sec. [n] are equal to one another. For /i 7^ to, we obtain 
non-isometric, distorted black holes. 

Another difficulty encountered in the numerical mod- 
eling of black holes is the wide discrepancy in length 



scales involved. The computational grid needs to be large 
enough to capture outgoing gravitational waves and to 
minimize boundary effects. The grid must also have high 
resolution in the interior to accurately resolve the strong 
gravitational fields of black holes. For a finite difference 
code, adaptive mesh refinement (AMR) is needed to sat- 
isfy both of these requirements: high resolution and a 
large grid. In this paper we introduce our AMR ellip- 
tic solver that allows us to solve accurately for puncture 
black hole data on very large grids. Evolution studies of 
these data sets are underway Q ■ 

In Sec. m we set up the equations defining the ini- 
tial value problem for a distorted puncture black hole. 
We show how the puncture data can be formulated to 
give isometric data sets, thus enabling a comparison with 
the results of Bernstein et al. In Sec. IIIII we give a 
brief description of our AMR elliptic solver. In Sec. IIVI 
we present sample results both for isometric and non- 
isometric black hole data. 



II. FORMULATION OF THE PROBLEM 

Following Bernstein et al. ((J 0, @ we write the line 
element for the physical metric <?y as 



ds 2 
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where q is a specified function of the spatial coordinates 
and tjj is the conformal factor. In this paper we restrict 
ourselves to the expression for q used in Ref . , namely, 
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where rj = ln(2r//i). In Eq. J5J), Qo, n, c and fi are 
adjustable constants that affect the size and type of dis- 
tortion. We also restrict ourselves to data sets defined at 
a moment of time symmetry. Thus, the extrinsic curva- 
ture vanishes and the momentum constraints are trivially 
satisfied. The Hamiltonian constraint reduces to 



V 2 i/> - -Ril) = 



(3) 



where V 2 and R are the Laplacian and scalar curvature 
of the conformal metric, defined by = ip~ 4 gij- 
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Solutions of Eq. on the manifold J? 3 arc Brill waves. 
The distorted black hole data sets of Bernstein et al. were 
obtained by solving Eq. (0) on the manifold Mx S 2 , with 
Robin conditions at the outer boundary [<9(r0— r)/dr = 
as r — ► oo] and isometry conditions at the inner boundary 
[d(^/rip)/dr = at r = fi/2]. We will obtain solutions of 
Eq. © on R x S 2 by following the puncture construction 
of Ref. 0. Thus, we write the conformal factor as ifi = 
u+m/(2r) and insert this into the Hamiltonian constraint 
© to obtain 



- 9 1 - m - 



s 



16r 



(4) 



This equation is solved for a continuous solution u on 
J? 3 with Robin boundary conditions d{ru — r)/dr = at 
r — > oo . Note that the "bare mass" m appears as a new 
parameter in the construction. 1 

We have had no difficulty in obtaining numerical so- 
lutions to Eq. l|Hl. However, from an analytical point 
of view, it is not immediately obvious to us whether or 
not Eq. (0} always admits solutions, and if so whether 
those solutions are unique. The problem of existence 
and uniqueness of solutions of Eq. on M 3 is similar to 
the problem of existence and uniqueness of solutions of 
Eq. on R 3 . There are two notable differences. First 
is the appearance of a "source" term mR/(\Qr) on the 
right-hand side of Eq. Q . Note that with the choice © 
for the function q, the scalar curvature R goes to zero 
rapidly, R ~ (i n ( r ))2 r in(A 1 /r)-2-2i„(2) at r o, so that 

R/r does not blow up at the puncture. The second differ- 
ence between equation (0J and the familiar initial value 
equation © is that, for we need not require the so- 
lutions to be positive. Rather, we would like to know 
if u is greater than — m/(2r) so that the combination 
ip = u + m/(2r) is everywhere positive. 

By construction the data sets of Bernstein et al. are 
isometric; that is, they are symmetric under reflections 
about the black hole throat r = fi/2. More precisely, 
each data set is described by a metric tensor whose com- 
ponents gij, as functions of r, 9 and </>, are unchanged by 
the coordinate transformation r — > r = /i 2 /(4r). We can 
solve for reflection symmetric data sets using the punc- 
ture method as well, simply by setting the parameters /i 
and m equal to one another: fj, = m. To see this, we first 
note that the Hamiltonian constraint on R x S 2 can be 
written as (see, for example, Appendix D of Ref. [Tol p 



V 2 - l -R) 
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where V 2 and R are the Laplacian and scalar curvature 
for the metric = c/ij/r 2 . Let "4>\{r) denote a solution 



1 The parameters fi and m are dimensionful. Thus, our data sets 
can be described in terms of the dimensionless parameters c, Qq, 
n, the dimensionless ratio n/m, and the dimensionless coordi- 
nates x/m, y/m, and z/m. 



of Eq. © , or equivalently Eq. (JSJ) , obtained by the punc- 
ture method. (For notational simplicity, we display only 
the r dependence in the solution ipi , and in ip2 below. In 
general these are functions of 9 and cj> as well.) This solu- 
tion ipi(r) has boundary behavior ipi( r ) ~ ► 1 as r ~~ * °°j 
and ip\ (r) — > m/(2r) as r — ► 0. Now observe that the line 
element ds 2 = gijdx 1 dx^ , with the function q chosen as in 
Eq. J2J), is invariant under reflections r — > f = /i 2 /(4r). 
The scalar operator V 2 — R/8 is invariant as well. By 
making the substitution r — > f = fi 2 /(4r) in Eq. © 
we see that ip2( r ), defined by y/rfair) = y/fipii^) with 
f = /i 2 /(4r), is also a solution of Eqs. JSJ) and ©. This 
solution satisfies the boundary conditions ^(r) — > m/fi 
as r — > oo and ^(r) — > l^/(2r) as r — > 0. If we choose 
/j = m, then the solutions V'i( r ) an d ^(f) satisfy the 
same equation and obey the same boundary conditions. 
If we assume that solutions to the puncture equation Q 
are unique, then the two solutions ipi( r ) an d ip^i 7 ") must 
in fact be identical. From the equality ipi(r) — ipiir) we 
find 



y/ripi (r) 



(6) 



r—fi 2 / (4r) 

The line element Q for this solution can be written as 

ds 2 = (^Mr)) 4 [e 2q (dr 2 /r 2 + d9 2 ) + sin 2 9 d<f> 2 ] . (7) 

Then the relation JSJl is the condition for the physical 
metric J7J to be invariant under the reflection r — > f = 
/V(4r). 

In Sec. IIVI we present results for the conformal factor 
tp for isometric (fi = m) and non-isometric (fi =/= m) 
data sets. For isometric data, we can check the reflection 
symmetry by comparing the ADM mass computed at the 
two infinities. At the "outer" infinity (r — > oo) the ADM 
mass is given by [TT1 | 



dS l V.;V 



(8) 



Numerically, the integral is computed at the outer bound- 
ary of the computational domain. The ADM mass at the 
"inner" infinity (r = 0) is found by expressing the metric 
JU in coordinates f, 9, tf>, where f = m 2 /(4r). We find 



ds 2 = 1 



^ 4 



2f / 



e 2q (df 2 +f 2 d9 2 ) + f 2 sin 2 9 d(/) 2 } , (9) 



where q and u are functions of r = m 2 /(4f), 9, and 4>. 
Provided q goes to zero sufficiently rapidly as f — ► oo, the 
metric is asymptotically flat at r = with ADM mass 
given by 



M = m u\ r=0 

The choice for 5 vanishes as r 
faster than any power of r. 



(10) 

like q - r ln ^/ r \ 
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III. NUMERICAL CODE 

Our numerical code solves Eq. using multigrid tech- 
niques with mesh refinement. We use the Paramcsh pack- 
age IT^ to implement mesh refinement and paral- 
lclization. Paramcsh divides the computational grid into 
blocks, with each block containing N zones. A block of 
data is refined by bisection — that is, the block is divided 
into eight "child" blocks (in three spatial dimensions) 
each containing N zones. 

Our code carries out multigrid V-cycles using the Full 
Approximation Storage (FAS) algorithm [l^ on non- 
uniform grid structures, with zone-centered data. It 
works with both fixed mesh refinement (FMR) and adap- 
tive mesh refinement (AMR). Working with FMR, we 
specify a non-uniform grid by hand. Working with AMR, 
we generally start with a coarse, uniform grid. The code 
V-cycles until the norm of the residual is less than the 
norm of the relative truncation error in each block. Any 
block whose relative truncation error is above a specified 
tolerance is flagged for refinement. Paramesh then re- 
builds the grid, the data is prolonged from the old grid to 
the new, and the V-cycle process begins again. Working 
in this mode takes the place of the Full Multigrid (nested 
iteration) Algorithm [l4j , where the solution on the final 
grid structure is reached by a succession of V-cycles of 
increasing peak resolution. 

We will provide details of the AMR-multigrid algo- 
rithm in a later publication, along with numerous code 
tests ITU. 



IV. RESULTS 

Bernstein et al. 0, 0, Q produced nonaxisymmetric 
distorted black holes using isometry conditions at the 
black hole throat. We have reproduced several of the 
initial data sets from Ref. using the puncture method, 
by setting /i = mas described in Sec. |H] As a check for 
isometry, we have confirmed that the ADM mass at the 
puncture, Mo, agrees with the ADM mass at infinity, 
Mqo, to several significant digits. 2 

Two tests were performed to confirm the consistency 
of our ADM mass calculations. For the first test we used 
a sequence of grids with fixed interior resolution but in- 
creasing outer boundary limits. We started by solving for 
the initial data on a grid with boundaries (in each dimen- 
sion) at —13 and 13, and having two additional box-in- 
box refinement levels ranging (in each dimension) from 
—6.5 to 6.5 and from —3.25 to 3.25. The resolution on the 
finest of the three levels was Ax = 0.05078. Other grids 



2 A direct numerical comparison of our ADM masses with those 
of Bernstein et al. is not possible, since their ADM masses are 
presented graphically. A visual comparison our ADM mass data 
with theirs shows good agreement. 
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FIG. 1: Contour plot of u in the z = plane for the case 
c = -2, Q = -0.5, n = 4, fi = m = 2. The full AMR grid 
extends from —100 to 100, while the inset shows the range 
-13 to 13. 



in the sequence were created by adding, one at a time, 
another fixed box-in-box refinement level while doubling 
the size of the computational domain. This test was in- 
tended primarily as a check on the sensitivity of the ADM 
mass at infinity to the location of the outer boundary. We 
found that the ADM mass M^ is unaffected to six digits 
and the ADM mass AI is unaffected to seven digits. For 
the second test we used a sequence of grids consisting of 
a fixed three-level box-in-box structure but increasing 
resolution throughout the computational domain. This 
test was intendend primarily to check the convergence 
properties of the ADM masses and Mo. The test 
showed that the ADM masses converge with second or- 
der accuracy. 

As a specific example let us consider the data set dis- 
cussed in Ref. ||| with c = —2, Qq = —0.5, n = 4, and 
fx = m = 2. The ADM masses for this data at infinity 
and at the puncture arc M^ = 2.2021 and M = 2.2027. 
Figure ^ shows a contour plot of the nonsingular part u 
of the conformal factor in the z = plane. The contour 
levels crossing the x axis between —100 and zero are, re- 
spectively, {1.0015,1.005,1.005,1.0015,0.96}. The con- 
tour levels crossing the y axis between —100 and zero are, 
respectively, {1.0015, 1.005, 1.01, 1.02, 1.1}. The inset in 
Fig-Dshows the range in x and y from —13 to 13. In the 
inset figure, the contours crossing the x axis between —13 
to zero are {1.005, 1.0015,0.96} and the contours cross- 
ing the y axis between —13 and zero are {1.01, 1.02, 1.1}. 
The function u forms two peaks along the y axis with 
maximum values 2.03, and two valleys along the x axis 
with minimum values 0.296. The value of u at the punc- 
ture is 1.10. The profiles of u along the three axes look 
similar to those shown in Fig. [5] 

With AMR, we are able to push the boundaries out 
quite far while maintaining high resolution around the 
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FIG. 2: The function u plotted along the three coordinate 
axes for a non-isometric distorted black hole with c = —2, 
Qo = —0.5, n = 4, /i = 2, m = 10. The thin solid line (top 
curve) is u along the j/-axis, the dotted line (middle curve) is 
u along the z-axis, and the thick solid line (bottom curve) is 
u along the z-axis. 

puncture. For the data described above, the computa- 
tional grid extends from —100 to 100 and has refined 
itself to reach a maximum relative truncation error of 
0.008. The lowest resolution region is equivalent to a 64 3 
grid with Ax = 3.125. The highest resolution region is 
equivalent to a 262, 144 3 grid with Ax = 0.000763. 

The puncture method allows us to define a new class of 
initial data: distorted black holes that are not isometric 
(/j 7^ m). Figure El shows such a data set with c = —2, 
Qo = —0.5, n = 4, n = 2, and m = 10. The results were 
computed on a grid extending from —100 to 100 in each 
dimension. The grid spacing for the lowest resolution re- 
gion, adjacent to the grid boundaries, is Ax = 3.125. The 
grid spacing for the highest resolution region, surround- 



ing the puncture, is Ax = 0.01221. The ADM mass at 
infinity is = 10.7, while the ADM mass at the punc- 
ture is Mq = 12.8. Figure El shows the behavior of the 
nonsingular part u of the conformal factor along the co- 
ordinate axes. The thin solid line plots u along the y 
axis, the thick solid line plots u along the x axis, and the 
dotted line plots u along the z axis. The two peaks in 
u along the y-axis have maximum values 4.26, and the 
two valleys along the x-axis have minimum values —1.38. 
The value of u at the puncture is 1.28. 

Note that for the non-isometric data described above, 
u takes on negative values in two small regions of radius 
~ 1 at locations ~ ±1 along the x axis. However, for this 
data set, the conformal factor ip = u + m/(2r) is every- 
where positive. We have explored other data sets that 
contain regions with u < 0, and in each case we find that 
the conformal factor ip is positive everywhere. For data 
with c = —2, Qo = —0.5, and n = 4, the most extreme 
data sets we have studied have ratios [ijvn = 1/300 and 
/i/m = 100. For the case /i/m = 1/300, the function 
u reaches a minimum value of ~ —120 but the confor- 
mal factor remains positive. For the case fi/m = 100 the 
function u, and therefore also ip, is always positive. 
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